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ABSTRACT 

This paper presents an improved k-e model and a second order closure model for low- 
Reynolds number turbulence near a wall. For the k-e model, a modified form of the eddy 
viscosity having correct asymptotic near- wall behavior is suggested, and a model for the pres- 
sure diffusion term in the turbulent kinetic energy equation is proposed. For the second order 
closure model, we modify the existing models for the Reynolds-stress equations to have proper 
near-wall behavior. A dissipation rate equation for the turbulent kinetic energy is also refor- 
mulated. The proposed models satisfy realizability and will not produce unphysical behavior. 
Fully developed channel flows are used for model testing. The calculations are compared with 
direct numerical simulations. It is shown that the present models, both the k-e model and the 
second order closure model, perform well in predicting the behavior of the near wall turbulence. 
Significant improvements over previous models are obtained. 


1. INTRODUCTION 

Many flow problems of practical importance such as diffusion controlled transition, drag- 
reduction, and flow separation and reattachment, are the result of low-Reynolds number tur- 
bulence near a wall. Accurate predictions of these flows require better understanding and 
modeling of turbulence near a wall. Most practical closure models (including k-e models and 
second order closure models) do not resolve a rather thin, viscosity-affected sublayer. They 
are usually incorporated with wall functions (Launder et a lJ 1 !, Shih and Lumleyl 2 !). However, 
wall functions are not suitable for many flows of practical importance where the equilibrium 
assumptions built in the wall functions are not valid. It is therefore necessary to extend and 
develop schemes for modeling turbulence directly down to the wall. There exist several low- 
Reynolds number k-e models which were discussed in detail by Patel et alJ 3 l. Patel et al. 
pointed out that the damping functions used in those models, especially the one for the eddy 
viscosity, need to be modified in order to improve model performance. Hanjalic and Launder^ 
proposed a Reynolds-stress closure model for near wall turbulence with no wall functions. How- 
ever, they did not use the full set of transport equations for the Reynolds-stress tensor in their 
implementation of the model. Most of the existing models that do not use wall functions do 
not properly predict the level of the normal components of the Reynolds-stress. For more 
discussions about the development of models for low-Reynolds number near wall turbulence, 
see Rodil 5 ! and Launder and Tselepidakis! 6 ! 

In this paper, the k-e model is presented in section 2. The near-wall asymptotic behavior 
of the eddy viscosity and the pressure transport term in the fc-equation is analyzed. According 
to that near-wall behavior, we proposed new models for the eddy viscosity and the pressure 
transport term. In addition, a modified dissipation rate equation is proposed following an 
argument similar to that of Lumleyl 7 !. In section 3, we present the second order closure model. 


* Work funded under Space Act Agreement C-99066-G. 


The near- wall behavior of turbulence (see e.g. Mansour, Kim & Moint 8 !, MKM hereafter) is used 
to form a set of model transport equations for the Reynolds-stress tensor and the dissipation 
rate of turbulent kinetic energy. We developed a near wall model for the combination of the 
velocity pressure-gradient correlation and dissipation rate tensor, and constructed a modeled 
dissipation rate equation. Near the wall, viscous effects become important and the reduction 
of velocity fluctuations normal to the wall becomes much more significant comparing to other 
velocity components. This wall effect makes the viscous diffusion terms (which are usually 
neglected in high-Reynolds number turbulence closures away from the wall) become one of 
the largest terms that must be properly balanced by the other terms in the Reynolds-stress 
equations. The proposed model equations satisfy the realizability condition which ensure no 
unphysical behavior. In section 4, we test the proposed models using a fully developed channel 
flow. For this flow, direct numerical simulation (Kim et al.I 9 J) data is available. This data was 
recently verified using three-dimensional particle tracking velocimeter measurements (Nishino 
and Kasagit 10 !). The modeled transport equations in this flow are one-dimensional and steady, 
hence model testing will be accurate. The results show that the k-e model proposed in this 
paper obtains an improvement over existing k-e model. The present second order closure 
model is capable of predicting wall behavior of various turbulent quantities. The mean velocity, 
shear stress and normal stresses calculated from the second order closure model are in a good 
agreement with the direct numerical simulation data. 


2. k-e MODEL 

2.1 Near- Wall Asymptotic Behavior 

To analyze the near-wall asymptotic behavior of the eddy viscosity and other turbulent 
quantities, we expand the fluctuating velocity and pressure in Taylor series about a normal 
distance y from the wall as follows: 

ui = b x y + Ci j/ 2 + diy 3 + . . . 

«2 = c 2 y 2 + d 2 y 3 + . . . 

, o , 3 (2.1) 

«3 = b 2 y + ay 1 + d 3 jr + . . . 

p = a p + b p y + c p y 2 + d p y 3 + ... 

where the coefficients a p ,6j t C2,... are functions of x,z and t. Based on the continuity and 
momentum equations, MKM show that the following relations between the coefficients hold: 

2 c 2 = — (&1 ,1 + 63,3) 

a Pt 1 = 2 i/ci 

a Pt 3 = 2 ^c 3 

b p - 2 vc 2 

where ( represents a derivative with respect to z,. The eddy viscosity is defined by 

-{u iUj ) = v T {Ui,j + Ujj) - | kiij (2.2) 

where ( ) stands for ensemble average and k = (u{U,)/2 = (q 2 )/ 2 is the turbulent kinetic 
energy. For plane shear flows, l>t = —(««)/ from Eq. (2.2). Using Eq. (2.1), we obtain the 
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(2.3) 


near-wall asymptotic behavior of the eddy viscosity: 

VT~zr~ — ~(b\C 2 )y 3 + ( — {^1 d ‘2 + C 1 C 2 ) + 2 {b\C 2 )(c\))y 4 + 0(y 5 ) 

oy 

That is, near the wall v t is 0(y 3 ), because dU/dy is usually 0(1). We shall see later that 
some existing models do not have this wall behavior. The near- wall asymptotic behavior of the 
turbulent kinetic energy k and its dissipation rate e = v{uijUij) are from Eq. (2.1) 

k = + ((b lCl ) + {b 3 c 3 ))y 3 + 0(y 4 ) (2.4) 

- = (b]) + (b 2 3 ) + 4({6 1 c 1 ) + {b 3 c 3 ))y + 0(y 2 ) (2.5) 

In the fc-equation, the pressure transport term II = becomes 

II = — 2u{(biCi) + ( b 3 c 3 ))y -f 0(y 2 ) (2.6) 

while the turbulent transport term, -( ku , ),<, is 0(y 3 ). Therefore, the pressure transport term 
becomes much more important than the turbulent transport term near the wall and must be 
properly modeled. 

2.2 Eddy Viscosity and Pressure Transport Term 

The eddy viscosity model can be in general written as vt = c u'C , where v! and i' are the 
turbulent characteristic velocity and length scale, respectively. In two-equation models, e.g. 
k-e models, k 1 / 2 is used as the characteristic velocity, and the length scale is characterized by 
ib 3 / 2 /f. Hence, the eddy viscosity is written as^ 13 ! 

J'T = C nf ft ~~ (2.7) 

where C ^ — 0.07, / M is a damping function. The form of the damping function is quite critical 
to the prediction of the mean flow field! 3 !. \ n fact, the prediction of the mean velocity depends 
on the eddy viscosity model. Therefore it is important for an eddy viscosity model to have 
proper near- wall behavior. We have examined the near- wall behavior of various eddy viscosity 
models which are listed in table 1. Table 1 shows that some of the k-e models do not have the 
correct near-wall behavior of the eddy viscosity. 

One expects that near the wall, the size of the large eddies (or energy containing eddies) 
should be order of the wall distance 0(y). Equations (2.4) and (2.5) show that fc 3 / 2 /e is 0(y 3 ). 
Hence, k 3 / 2 /e is not an appropriate quantity to represent the length scale of the large eddies 
near the wall. However, we can define a variable e = e — i/k f ik t {/2k which has a nice property: 
e approaches e away from the wall and is 0(y 2 ) near the wall. Therefore k 3 ! 2 fe is 0(y ) and is 
a proper quantity to characterize the length scale of the large eddies. With this length scale, 
the eddy viscosity model should be written as 

= Cpfp— ( 2 - 8 ) 

e 

Now, in order for v T to have correct near- wall behavior, the damping function must be O(y) 
near the wall and approach 1 away from the wall. The damping functions used in various k-€ 
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models are also listed in table 1. If we consider the presence of the wall as the main effect on 
the eddy viscosity, then we can assume is mainly a function of (defined as u T y/v , where 

u T is the friction velocity). The form of can be determined quite accurately if we know k 

and e, for example, from the direct numerical simulations. We may also optimize the following 
simple form by numerical experiments: 

fn = 1 - exp(— a x y + - a 2 y +2 - a 3 y +3 - a 4 y + 4 ) (2.9) 

The optimal values for channel flows are a x = 6 x 10 _3 ,a 2 = 4 x 10~ 4 ,a3 = -2.5 x 10 -6 ,a 4 = 
4 X lO -9 . They can be further tuned using the direct numerical simulation data. 

For the term of the pressure transport of the turbulent kinetic energy, we propose: 

tt r 0.05 vt , x 

1/ M [1 -exp(-j/+)] a k ’Ti ( 2 - 10 ) 

which has a similar form to the standard turbulent transport model, but with a coefficient to 
ensure the near-wall behavior, Eq. (2.6). Here the coefficient 0.05 is an optimal model constant 
for the channel flows. 


2,3 Modeled k-c Equation 


To complete the eddy viscosity model, we need the modeled equations for the turbulent 
kinetic energy and for its dissipation rate. In general one may write the modeled A; -equation 
as follows: 


dk dk 

Hi + Uj ~ = 


dx , 


d 

dxj 


r f "T.dk dUi OUj.dUi 

K*' + — )— ] + n + u t{— + + D 


dxj 


dxj dxi uxj 


(2.11) 


where D is zero in the present k-t model, because c is finite at the wall. In k-e models that 
assume e = 0 at the wall, an artificial term D must be added in order to balance the viscous 
diffusion term at the wall. The form of D for various k-€ models is listed in table 1. 


For the modeled dissipation rate equation, we write 


(,t + Uj€j = [(1 -(- 1 'T/<r ( )c,j\,j - (2.12) 

where ’J 1 stands for the entire mechanism of the production and destruction of the dissipation 
rate e. At the level of the k-e model, we assume that $ is a function of u, u j, k, e, e, Uij 
and Uijk- Because $ is an invariant, it must be a function of the invariants that can be 
constructed from at least following quantities: R t , i/ T Ui t jUij/e and v^TU itjk U i}jk k/ee, where 
R t is the turbulent Reynolds number k 2 /ve. We now expand $ in a Taylor series about these 
invariants and take only the linear terms. We obtain 


^ + i>i 1 - 'l J + favvrUijkUi, ,j k — (2.13) 

e ee 

where the coefficients ip o> V’l an ^ 1 P 2 are in general functions of R t . Finally, we write the 
dissipation rate equation as 


de de 

di + Uj d 71 


+ Tp-tp + - C.A f + £ 


(2.14) 
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where C\ and C 2 are model constants, and f\ and fi are functions of R t . The term E in the 
present model comes from the last term in Eq. (2.13): E = vi/xUijkUijk, where we have taken 
= —1. The forms of E and C*i, C 2 , f\ and fi for various k-e models are listed in table 1. 


3. SECOND ORDER CLOSURE MODEL 
3.1 Reynolds-Stress Equation 

For an incompressible flow, the Reynolds-stress equations can be written 

u i u j ) “ Pij + Tij + E\j + II ij — Sij (3.1) 

where ( ) stands for an ensemble average, D/Dt = d/dt+ Ukd/dxk, and the terms on the right 
hand side of the above equations are identified as follows: Pij = — ( UiUk)Uj t k — 

Pij = (Uj'Ujf Ufc) , E ij — v(U{Uj ) , II * j — “(^ tP t j T ^ jP ,i) 1 £ij = {^i y k^j ,k) ' 1^ 

Eq. (3.1), the turbulent transport Ta, the velocity pressure-gradient correlation II ij and the 
dissipation rate Sij must be modeled. Near-wall analysis of existing models shows that they do 
not exhibit a proper near- wall behavior, which indicates that the modeled transport equations 
are not properly balanced near the wall. This may explain the poor quantitative predictions 
from existing models for near wall turbulence, specifically, in predicting the normal stresses d 6 3 

To analyze the near- wall asymptotic behavior of the different terms in the Reynolds-stress 
equations, we use (2.1) to estimate each term in the Reynolds-stress equations. For example, 
for the (^1^2) component, 


D{u\U 2 ) 

Dt 

P12 

Tn 

E\2 

IIl2 


= 0 ( 2 / 3 ) 

= O0/ 4 ) 

= OG/ 4 ) 

= 61/(61 c 2 )y + 0(y 2 ) 

= -2u{b x c 2 )y + 0{y 2 ) 


e 12 = 4i/(bic 2 )y + 0(y 2 ) 


(3.2) 


The near-wall budget for other components can be expressed similarly. We see from these 
budgets that the viscous diffusion term Dij is always the largest term. Its leading term is 
balanced by the dissipation rate £{j or by the combination of the dissipation rate and the 
velocity pressure-gradient correlation, -E{j + n,-j. In this paper, we derive a near wall model 
for £ij and II y in a rational way, using the wall behavior as a model constraint. Away from 
the wall, we adopt the high-Reynolds number realizable model of Shih and LumleyJ 11 ! Near 
the wall, we write the combination velocity pressure-gradient and dissipation rate as 


£ ij - 


(3-3) 


Note, $ij is not trace free because it includes the dissipation rate and the pressure diffusion. 
We model <£ t j by assuming that it is a linear function of {u{Uj), and function of a unit vector 
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normal to the surface, n,*. Under these assumptions, the most general form of which satisfies 
symmetry is 

$ij = 0‘00‘^i'fij + cLo^ij + a\(uiUj) + a,2({uiUk)njrik + {ujU^niTi^) 

+ a 3 {ukUi)n k niniTij + a 4 (u k ui)n k ni6ij 

Using Eq. (3.4) and the wall constraint, 

■^0 "h ^ 2 y^«j = 0 as y — ► 0 (3.5) 

the coefficients can be determined as a 00 = a 0 = 0, O! = 2, a 2 = 4, and a 3 = 2, with an 
undetermined a 4 , which is set to be zero according to numerical experiments. The above model, 
Eq. (3.4), is valid only in the region near the wall. The influence of the wall should vanish at 
the place far away from the wall. Numerical experiments suggest the use of a damping function 
of the form f w = exp {-(R t /C) 2 ), where R t = ( 9 2 ) 2 /9i/e, C = 1.358 R° er 4 \ R eT = u r 6/u, u T is 
the friction velocity, and 6 is the thickness of the boundary layer or half width of a channel. 
Therefore, the model in Eq. (3.3) can be written as follows: 

flij — Eij = — f w j^[2(u{Uj) + 4({uiu k )rijn k + ( UjU k )nin k ) + 2{u k ui)n k n[Tiinj] (3.6) 


Awa y fr° m the wall > the velocity pressure-gradient correlation is split into a rapid part 
II- and a slow part IljA Lumley’s high-Reynolds number model! 7 ! for the return to isotropy 

term II— ~ Eij has been successfully used in conjunction with wall functions. However, this 
return to isotropy model does not vanish as y — ► 0, and should be modified with a damping 
function if one wants to use it to the wall. We use the same f w damping function to minimize 
the arbitrary choice of the damping functions, and set 

n <? - + f «*;)(1 - /«). (3-7) 

where b{j = ( UiUj)/(q 2 ) - 8ij/ 3. The rapid part of the velocity pressure-gradient, II A , is 
modeled using the model of Shih and Lumleyt 2 ! with a modified coefficient C. Notice that 
this rapid model decays as 0(y 2 ) as y — » 0. Therefore, we use it down to the wall without a 
damping function. Shih and Lumley!^! set C ~ 0.8 when the model was used in conjunction 
with wall functions. Here, we modify the coefficient as a function of the turbulent Reynolds 
number: C = 0.8[1 - exp(-(f? t /40) 2 )]; C will approach 0.8 when R t becomes large. 

For the third moments, we use the model of Hanjalic and LaunderW, but with a different 
model constant: 


Iql \ 

(uiUjUk) = -0.07— •«,■), p + (u J -u J ,)(tt j u*) i|l + (uiU p )(ujU k )' P } (3.8) 
3.2 Dissipation Rate Equation 

We write the dissipation rate equation as 


c,t + Ui€,i - = 


(e 

¥) 




(3.9) 
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At the level of the second order closure, we assume that $ is a function of (uiuj), c, v, Uij 
and Uijk- Since ^ is an invariant, it should be at least a function of //, III , R t , ,J , 

UiJkUiJl{ u k u l)/ UiJkUj,il{ U k u l)/ Jq?j i an( ^ U\,jkU l,ij{ u k u {) / Jq2) ^ W ^ ere 

II and III are the invariants of 6,j. We then expand $ in Taylor series about the invariants 
and take only first order linear terms. We obtain 


¥ = -00 + 01 


(ujU j ) Ujj 

e 


, UijkUiji('U'kV'i) 

‘ g g _g 

(<? 2 > (<7 2 ) ^ 


^ 6 6 _€ T g € g. 

<? 2 > <<? 2 ) * <? 2 > (F) 1 / 


(3.10) 


where the coefficients 0o, 0i, 02 1 03 an d 04 are, in general, functions of 77, /// and R t . Now 
if we consider the terms in (3.10) containing the second derivative of the mean velocity as a 
model of the term —2 v(ujUi tk )Ui }k j (the term in the exact equation of e containing the second 
derivative of the mean velocity) and use continuity, {uju kik ) = 0, we obtain 03 = 0, 04 = — 02. 
In addition, at high Reynolds number or when the small scales are isotropic, the dissipation 
rate should not respond instantaneously to a change in the mean velocity derivatives. Under 
this constraint, we require 02 to vanish as the turbulence becomes isotropic. The simplest 
form that satisfies this condition is 09 = —0.15(1 — 7), where Lumley’s parameter* 7 ! F = 
1 + 9b{jbjkbki — 4.5 bijbij approaches one as the turbulence becomes isotropic. Finally, the 
modeled dissipation rate equation can be written as 


C,t + Uit yi - - (eui)) t i - ip 0 


€€ 

¥) 


£ v(q 2 ) 

- Vij-2j(uiUj)Uij - rp2—j~(ukUi)(Uiji - Ui tij )Uijk 


(3.11) 


where ipo, and ip\ are Lumley’s coefficients.^ To close the above equation, we use Hanjalic and 
Launder’s diffusion model with a different model constant: 


/ (F \ 

(eu k ) = “0*07 : ~{u k u p )e tP 


(3.12) 


4. FULLY DEVELOPED TURBULENT CHANNEL FLOW 

To test the models developed in this paper, we chose a fully developed channel flow. The 
Reynolds number based on the friction velocity and the half width of the channel Re r is 180, 
for which direct numerical simulation* 9 ! and experimental data* 10 ! are available for comparison. 
The Reynolds-stress equations for this flow are one-dimensional and steady, therefore the model 
testing will be accurate. 

The main results of the k-e models are shown in figures (1) to (4). These figures present 
predictions for the mean velocity, turbulent kinetic energy and its dissipation rate and eddy 
viscosity profiles, using several k-e models* 13 !’* 14 !’* 15 !’* 16 ! including the present one. These re- 
sults show clearly that the present k-e model gives significant improvements over previous k-e 
models according to the comparisons with the direct numerical simulations. The main results 
of the present second order closure model are shown in the figures (5) to (8). All the predic- 
tions are in a good agreement with the direct numerical simulation data. Figures (5) and (6) 
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present the comparisons of predicted mean velocity and turbulent kinetic energy using different 
models (both k-c models and second order closure models) with the numerical simulation data. 
The comparison with data of the normal components of the stresses are shown in figure (7). 
Predictions for all the components are quantitatively accurate. The near wall behavior of the 
three normal components is also reflected by the present models. None of the previous models 
have succeeded in predicting these trends.l 5 )’! 6 ! 

It should be noted, in our comparison with the data, the channel flow was computed 
at fairly low Reynolds number. At high Reynolds numbers, the models still predict the log 
law profile for the mean velocity. Comparison with new channel data! 18 ! (see figure (8)) at 
Re T = 393 shows that the normal stresses near the wall are still well predicted. However, in 
the core region the model slightly overpredicts the level of the turbulent kinetic energy. 
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TABLE L Constants and functions for various K — e models 


Model 

c w -BC 

C, 

c tl 

C« 2 

Ok 

o f . 

JL 

0 

.09 

1.45 

2.0 

1.0 

1.3 

Reynolds 

u^K 

V dy 3 

.084 

1.0 

1.83 

1.69 

1.3 

LB 

i/SlK 

V dy 2 

.09 

1.44 

1.92 

1.0 

1.3 

Chien 

0 

.09 

1.35 

1.8 

1.0 

1.3 

NH 

0 

.09 

1.45 

1.9 

1.0 

1.3 

Present 

dy 3 

.09 

1.45 

2.0 

1.3 

1.3 


Model 


/ 1 

/ 2 

JL 

“PtTO) 

1.0 

1 - .3exp(-i? 2 ) 

Reynolds 

1 - exp(-.0198#*:) 

1.0 

[l-.3exp (-Rj/9))/^ 

LB 

[1 - exp(-.0165/Zjt)] 2 

x(l + ¥) 

l + (f) 3 

1 - exp (~R 2 t ) 

Chien 

1 - exp(-.0115y + ) 

1.0 

1 - .22 exp( — / 36) 

NH 

[1 - exp(-i/ + /26.5)] 2 

1.0 

1 - .3exp(-i2 2 ) 

Present 

Eq. (2.9) 

1.0 

1 - .22 exp( — i? 3 /36) 


Model 

n 

D 

JL 

0 


Reynolds 

0 

0 

LB 

0 

0 

Chien 

0 

2uF< 

y 2 

NH 

0 

-H^) 2 

Present 

Eq. (2.10) 

0 


Rt = K 2 jvc, Rk - y/Ky/u, = u r yju. 


E i/x 

2 -t(0) 2 0(</ 3 ) 

0 0(j/ 5 ) 


o 0(j/ 4 ) 

— exp( — .5j/ + ) 0(y 3 ) 

^t(0) 2 0(2/ 3 ) 
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Figure 1. Comparisons between k-e models and 
DNSM-data for mean velocity profiles in the 2D channel. 
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Figure 2. Comparisons between fc-t models and 
DNSl°l-data for turbulent kinetic energy in the 2D channel. 
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Figure 3. Comparisons between k-e models and 
DNSl°l-data for the dissipation rate in the 2D channel. 
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Figure 4. Comparisons between k-e models and 
DNSl°l-data for eddy viscosity in the 2D channel. 
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Figure 6. Turbulent kinetic energy profiles, Re r = 180. 
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Figure 7. RMS of fluctuating velocities, Re r — 180. 

Lines represent present model; Symbols are DNSl 0 l-data. 
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Figure 8. RMS of fluctuating velocities, iie r =393. 

Lines represent present model; Symbols are DNSl 16 l-data. 
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